import numpy as np
from EF2D.InpReader import read_abaqus_inp
from EF2D.mesh import MeshBlock
from EF2D.Elements import PS4
from EF2D.Loads import *
from EF2D.node import Node
from EF2D.fea import FEA
from EF2D.tools import *
from icecream import ic
from typing import List,Dict
from PostPlot.pyVersion.VTKplot import meshplot

initNumpyPrintOptions()
# define material properties[1100铝合金 ANSYS LS-DYNA UserGuide ｀ANSYS'2008.]
E=69000 #MPa
nu=0.3
t=1
pho=2.71e-9 #t/mm^3
Xforce=1000 #N

def test_PatchTest_PS4():
    # 读取网格数据
    nodes, elements, nsets, _ = read_abaqus_inp(".\\tests\\data\\20250228PatchTestPS4.inp")
    ic(nodes)
    ic(elements)
    ic(nsets)
    
    # create meshblock
    NodesList :Dict[int,Node] = {int(n[0]):Node(n[0],x=n[1],y=n[2],z=n[3]) for n in nodes}
    ElemsList:List[PS4]=[PS4(label=e[0],
                    node1=NodesList[e[1]],
                    node2=NodesList[e[2]],
                    node3=NodesList[e[3]],
                    node4=NodesList[e[4]],
                    E=E,nu=nu,thickness=t) for e in elements]
    MbPS4=MeshBlock(name="PS4",elements=ElemsList)
    
    # create loads
    LoadsList:List[ConcentratedLoad]=[ConcentratedLoad("Xforce",nodeLabels=nsets['loadx'],dofth=1,value=Xforce),]
    # create bc,fix Left Nodes
    BcList:List[DisplacementBc]=[DisplacementBc("fixed",nodeLabels=nsets['fix'].astype(int),dofth=1,value=0.0),
                                    DisplacementBc("fixed",nodeLabels=nsets['fix'].astype(int),dofth=2,value=0.0),]
    # create FEA object
    fea=FEA(meshblocks=[MbPS4,],loads=LoadsList,bcs=BcList)
    # solve
    # fea.showKg()
    # fea.showMg()
    
    # 选择求解类型为 静力分析
    fea.Solve(StepType='Static')
    ic(fea.get_U().T)
    ic(fea.get_F_react().T)
    
    # ploting
    meshplot(nodes=nodes[:,1:],
                elements=[{'NodeIndices':e[1:].tolist(),'CellType':PS4.VTK_TYPE} for e in elements],
                scalars=fea.get_U()[1::2,0],
                range=None,
                show_lines=True,
                NumColors=10,
                title='Mesh with Uy')
    meshplot(nodes=nodes[:,1:],
                elements=[{'NodeIndices':e[1:].tolist(),'CellType':PS4.VTK_TYPE} for e in elements],
                scalars=fea.get_U()[0::2,0],
                range=None,
                show_lines=True,
                NumColors=10,
                title='Mesh with Ux')
    
if __name__ == '__main__':
    # 位移结果与ABAQUS 基本一致. abaqus 结果为:2.84e-2, 程序结果为:2.82e-2
    test_PatchTest_PS4()